%%%%% Import Estimates
 
if Net_OUT==0
    load(['PE_Estimates_',num2str(sim_tol),'_OR.mat'])
    load(['Estimates_NoA_M1_OR.mat'])
    Chain_Out_M1 = Chain_Out_NoA;
    load(['Estimates_NoA_M2_OR.mat'])
    Chain_Out_M2 = Chain_Out_NoA;    
else
    load(['PE_Estimates_',num2str(sim_tol),'_OUT.mat'])
    load(['Estimates_NoA_M1_OUT.mat'])
    Chain_Out_M1 = Chain_Out_NoA;
    load(['Estimates_NoA_M2_OUT.mat'])
    Chain_Out_M2 = Chain_Out_NoA;       
end
 
 
if T2C == 0
    load(['Mats_T1.mat'])
    Cons = Cons_T1;
    NetOut = NetOut_T1;
else
    load(['Mats_T2C.mat'])
    Cons = Cons_T2C;
    NetOut = NetOut_T2C;    
end
 
obs = size(NetOut.G_ij,1);
 
 
BGD_hat_M1 = mean(Chain_Out_M1.BGD_hat_b,2);
V_C_M1 = mean(Chain_Out_M1.V_C_b,2);
V_C_M1 = [V_C_M1, [V_C_M1(2); V_C_M1(1)]];
Est_M1 = struct('BGD_hat',[BGD_hat_M1(1:7); zeros(7,1); 0; BGD_hat_M1(8:13); zeros(34,1)],'V_C',V_C_M1,'A_hat',zeros(Cons.n,1), ...
    'l_M_hat',zeros(Cons.n,1),'beta_M_hat',[mean(Chain_Out_M1.beta_M_hat_b,2);zeros(7,1)],'sigma2_M_hat',mean(Chain_Out_M1.sigma2_M_hat_b), ...
    'alpha_educ_m1',mean(Chain_Out_M1.alpha_educ_m1_b,2),'alpha_educ_m3',mean(Chain_Out_M1.alpha_educ_m3_b,2), ...
    'alpha_gender_m1',mean(Chain_Out_M1.alpha_gender_m1_b,2),'alpha_gender_m3',mean(Chain_Out_M1.alpha_gender_m3_b,2), ...
    'Sigma2_educ_m1',mean(Chain_Out_M1.Sigma2_educ_m1_b,2),'Sigma2_educ_m3',mean(Chain_Out_M1.Sigma2_educ_m3_b,2), ...
    'Sigma2_gender_m1',mean(Chain_Out_M1.Sigma2_gender_m1_b,2),'Sigma2_gender_m3',mean(Chain_Out_M1.Sigma2_gender_m3_b,2));
 
 

BGD_hat_M2 = mean(Chain_Out_M2.BGD_hat_b,2);
V_C_M2 = mean(Chain_Out_M2.V_C_b,2);
V_C_M2 = [V_C_M2, [V_C_M2(2); V_C_M2(1)]];
Est_M2 = struct('BGD_hat',[BGD_hat_M2(1:14); 0; BGD_hat_M2(15:27); zeros(27,1)],'V_C',V_C_M2,'A_hat',zeros(Cons.n,1), ...
    'l_M_hat',zeros(Cons.n,1),'beta_M_hat',mean(Chain_Out_M2.beta_M_hat_b,2),'sigma2_M_hat',mean(Chain_Out_M2.sigma2_M_hat_b), ...
    'alpha_educ_m1',mean(Chain_Out_M2.alpha_educ_m1_b,2),'alpha_educ_m3',mean(Chain_Out_M2.alpha_educ_m3_b,2), ...
    'alpha_gender_m1',mean(Chain_Out_M2.alpha_gender_m1_b,2),'alpha_gender_m3',mean(Chain_Out_M2.alpha_gender_m3_b,2), ...
    'Sigma2_educ_m1',mean(Chain_Out_M2.Sigma2_educ_m1_b,2),'Sigma2_educ_m3',mean(Chain_Out_M2.Sigma2_educ_m3_b,2), ...
    'Sigma2_gender_m1',mean(Chain_Out_M2.Sigma2_gender_m1_b,2),'Sigma2_gender_m3',mean(Chain_Out_M2.Sigma2_gender_m3_b,2));



 
V_C_WithA = mean(PE_Chain_Out.V_C_b,2);
V_C_WithA = [V_C_WithA, [V_C_WithA(2); V_C_WithA(1)]];
 
Est_WithA = struct('BGD_hat',mean(PE_Chain_Out.BGD_hat_b,2),'V_C',V_C_WithA,'A_hat',zeros(Cons.n,1), ...
    'l_M_hat',zeros(Cons.n,1),'beta_M_hat',mean(PE_Chain_Out.beta_M_hat_b,2),'sigma2_M_hat',mean(PE_Chain_Out.sigma2_M_hat_b), ...
    'alpha_educ_m1',mean(PE_Chain_Out.alpha_educ_m1_b,2),'alpha_educ_m3',mean(PE_Chain_Out.alpha_educ_m3_b,2), ...
    'alpha_gender_m1',mean(PE_Chain_Out.alpha_gender_m1_b,2),'alpha_gender_m3',mean(PE_Chain_Out.alpha_gender_m3_b,2), ...
    'Sigma2_educ_m1',mean(PE_Chain_Out.Sigma2_educ_m1_b,2),'Sigma2_educ_m3',mean(PE_Chain_Out.Sigma2_educ_m3_b,2), ...
    'Sigma2_gender_m1',mean(PE_Chain_Out.Sigma2_gender_m1_b,2),'Sigma2_gender_m3',mean(PE_Chain_Out.Sigma2_gender_m3_b,2), ...
    'alpha_educ_m2',mean(PE_Chain_Out.alpha_educ_m2_b,2),'alpha_educ_m4',mean(PE_Chain_Out.alpha_educ_m4_b,2), ...
    'alpha_gender_m2',mean(PE_Chain_Out.alpha_gender_m2_b,2),'alpha_gender_m4',mean(PE_Chain_Out.alpha_gender_m4_b,2), ...
    'Sigma2_educ_m2',mean(PE_Chain_Out.Sigma2_educ_m2_b,2),'Sigma2_educ_m4',mean(PE_Chain_Out.Sigma2_educ_m4_b,2), ...
    'Sigma2_gender_m2',mean(PE_Chain_Out.Sigma2_gender_m2_b,2),'Sigma2_gender_m4',mean(PE_Chain_Out.Sigma2_gender_m4_b,2));
 
Sims_Out = struct('y1_educ_M1_m3',[],'y1_gender_M1_m3',[],'y1_educ_M2_m3',[],'y1_gender_M2_m3',[], ...
    'y1_educ_m3',[],'y1_gender_m3',[],'y1_educ_m4',[],'y1_gender_m4',[], ...
    'P_sim',[],'P_random',P_random,'sim_tol',sim_tol);
 
Cons.NETMISS = ones(obs,1);
Cons.A_MISS = ones(Cons.n,1);
Cons.y1_educ_MISS = ones(Cons.n,1);
Cons.y1_gender_MISS = ones(Cons.n,1);
 
y1_start = struct('y1_educ_m1',zeros(Cons.n,1),'y1_educ_m2',zeros(Cons.n,1),'y1_educ_m3',zeros(Cons.n,1),'y1_educ_m4',zeros(Cons.n,1), ...
    'y1_gender_m1',zeros(Cons.n,1),'y1_gender_m2',zeros(Cons.n,1),'y1_gender_m3',zeros(Cons.n,1),'y1_gender_m4',zeros(Cons.n,1));
 
 
for i=1:Sim_reps
    
    i
    Cons_rep = Cons;
    %%%%% When called for, randomly assign P
    if P_random == 1
        % Impute P
        J = rand(Cons_rep.n,1);
        if L_educ_pref==1
            J = J + (ones(Cons_rep.n,1)-Cons_rep.L_educ);
        end
        if M_educ_pref==1
            J = J + (ones(Cons_rep.n,1)-Cons_rep.M_educ);
        end        
        if H_educ_pref==1
            J = J + (ones(Cons_rep.n,1)-Cons_rep.H_educ);
        end
        if L_gender_pref==1
            J = J + (ones(Cons_rep.n,1)-Cons_rep.L_gender);
        end
        if M_gender_pref==1
            J = J + (ones(Cons_rep.n,1)-Cons_rep.M_gender);
        end        
        if H_gender_pref==1
            J = J + (ones(Cons_rep.n,1)-Cons_rep.H_gender);
        end
 
        
        P_rep = [];
        z0 = 0;
        for s=1:Cons_rep.schools
            J_s = J(z0+1:z0+Cons_rep.sch_size(s),:);
            [J_1, J_2, J_3] = unique(J_s,'sorted');
            P_rep_s = max(ceil((13.5-J_3)/13),sparse(Cons_rep.sch_size(s),1));
            P_rep = [P_rep; P_rep_s];
            z0 = z0 + Cons_rep.sch_size(s);
       end
      
        
        % Reset data matrics
        X_i_rep = Cons_rep.I_i*[Cons_rep.X_raw, P_rep, bsxfun(@times,P_rep,Cons_rep.X_raw)];
        X_j_rep = Cons_rep.I_j*[Cons_rep.X_raw, P_rep, bsxfun(@times,P_rep,Cons_rep.X_raw)];
        X_j_doti_rep = Cons_rep.Trans_doti*X_j_rep;
 
        mean_Xk_less_ij_rep = bsxfun(@rdivide,Cons_rep.I_i*Cons_rep.I_i'*X_j_rep - X_j_rep,Cons_rep.I_i*Cons_rep.I_sum');
        X_j_mean_Xk_less_ij_rep = X_j_rep.*mean_Xk_less_ij_rep;
        X_i_mean_Xk_less_ij_rep = X_i_rep.*mean_Xk_less_ij_rep;
 
        Cons_rep.X_i = X_i_rep;
        Cons_rep.X_j = X_j_rep;
        Cons_rep.P = P_rep;
        Cons_rep.X_j_doti = X_j_doti_rep;
        Cons_rep.mean_Xk_less_ij = mean_Xk_less_ij_rep;
        Cons_rep.X_j_mean_Xk_less_ij = X_j_mean_Xk_less_ij_rep;
        Cons_rep.X_i_mean_Xk_less_ij = X_i_mean_Xk_less_ij_rep;
    else
        P_rep = sparse(Cons_rep.P);
    end
 
    
    %%%%% Model 1: Simple (without A's, without P's)
    % (1) Impute network
    G_ij_start = zeros(size(NetOut.G_ij,1),1);
    [G_ij_M1_rep, A_NoA_rep] = Fn_ImputeNetwork_Mar2022(G_ij_start, Est_M1, Cons_rep, sim_tol, 0);
 
    if Net_OUT==0
        W_ij_M1 = exp(G_ij_M1_rep) + exp(Cons_rep.Flip_ij*G_ij_M1_rep);
    else
        W_ij_M1 = exp(G_ij_M1_rep);
    end
    W_ij_M1 = Cons_rep.W_trans*W_ij_M1;
    W_ij_M1_rep = zeros(Cons_rep.n,Cons_rep.n);
    k = 0;
    for s=1:Cons_rep.schools
        W_ij_M1_rep(k+1:k+Cons_rep.sch_size(s),k+1:k+Cons_rep.sch_size(s)) = reshape(W_ij_M1(k+1:k+Cons_rep.sch_size(s)^2,1),Cons_rep.sch_size(s),Cons_rep.sch_size(s));
        k = k + Cons_rep.sch_size(s);
    end
    W_ij_imp_M1 = bsxfun(@rdivide,W_ij_M1_rep,sum(W_ij_M1_rep,2));
 
 
    % (2) Impute Outcomes
    [y1_NoA_imp_M1] = Fn_ImputePE_Mar2022(Est_M1, zeros(Cons_rep.n,1), W_ij_imp_M1, Cons_rep, y1_start);
    
    Sims_Out.y1_educ_M1_m3 = [Sims_Out.y1_educ_M1_m3, y1_NoA_imp_M1.y1_educ_m3];
    Sims_Out.y1_gender_M1_m3 = [Sims_Out.y1_gender_M1_m3, y1_NoA_imp_M1.y1_gender_m3];        
    
    
            
    %%%%% Model 2: Simple (without A's)
    % (1) Impute network
    G_ij_start = zeros(size(NetOut.G_ij,1),1);
    [G_ij_M2_rep, A_NoA_rep] = Fn_ImputeNetwork_Mar2022(G_ij_start, Est_M2, Cons_rep, sim_tol, 0);
 
    if Net_OUT==0
        W_ij_M2 = exp(G_ij_M2_rep) + exp(Cons_rep.Flip_ij*G_ij_M2_rep);
    else
        W_ij_M2 = exp(G_ij_M2_rep);
    end
    W_ij_M2 = Cons_rep.W_trans*W_ij_M2;
    W_ij_M2_rep = zeros(Cons_rep.n,Cons_rep.n);
    k = 0;
    for s=1:Cons_rep.schools
        W_ij_M2_rep(k+1:k+Cons_rep.sch_size(s),k+1:k+Cons_rep.sch_size(s)) = reshape(W_ij_M2(k+1:k+Cons_rep.sch_size(s)^2,1),Cons_rep.sch_size(s),Cons_rep.sch_size(s));
        k = k + Cons_rep.sch_size(s);
    end
    W_ij_imp_M2 = bsxfun(@rdivide,W_ij_M2_rep,sum(W_ij_M2_rep,2));
 
 
    % (2) Impute Outcomes
    [y1_NoA_imp_M2] = Fn_ImputePE_Mar2022(Est_M2, zeros(Cons_rep.n,1), W_ij_imp_M2, Cons_rep, y1_start);
    
    Sims_Out.y1_educ_M2_m3 = [Sims_Out.y1_educ_M2_m3, y1_NoA_imp_M2.y1_educ_m3];
    Sims_Out.y1_gender_M2_m3 = [Sims_Out.y1_gender_M2_m3, y1_NoA_imp_M2.y1_gender_m3];        
                          
    
    %%%%% Models 3 and 4: With A's
    % (1) Impute network
    G_ij_start = zeros(size(NetOut.G_ij,1),1);
    [G_ij_WithA_rep, A_WithA_rep] = Fn_ImputeNetwork_Mar2022(G_ij_start, Est_WithA, Cons_rep, sim_tol, 0);
 
    if Net_OUT==0
        W_ij_WithA = exp(G_ij_WithA_rep) + exp(Cons_rep.Flip_ij*G_ij_WithA_rep);
    else
        W_ij_WithA = exp(G_ij_WithA_rep);
    end
    W_ij_WithA = Cons_rep.W_trans*W_ij_WithA;
    W_ij_WithA_rep = zeros(Cons_rep.n,Cons_rep.n);
    k = 0;
    for s=1:Cons_rep.schools
        W_ij_WithA_rep(k+1:k+Cons_rep.sch_size(s),k+1:k+Cons_rep.sch_size(s)) = reshape(W_ij_WithA(k+1:k+Cons_rep.sch_size(s)^2,1),Cons_rep.sch_size(s),Cons_rep.sch_size(s));
        k = k + Cons_rep.sch_size(s);
    end
    W_ij_imp_WithA = bsxfun(@rdivide,W_ij_WithA_rep,sum(W_ij_WithA_rep,2));
 
    % (2) Impute Outcomes
    [y1_WithA_imp] = Fn_ImputePE_Mar2022(Est_WithA, A_WithA_rep, W_ij_imp_WithA, Cons_rep, y1_start);   
    
    Sims_Out.y1_educ_m3 = [Sims_Out.y1_educ_m3, y1_WithA_imp.y1_educ_m3];
    Sims_Out.y1_educ_m4 = [Sims_Out.y1_educ_m4, y1_WithA_imp.y1_educ_m4];
    
    Sims_Out.y1_gender_m3 = [Sims_Out.y1_gender_m3, y1_WithA_imp.y1_gender_m3];
    Sims_Out.y1_gender_m4 = [Sims_Out.y1_gender_m4, y1_WithA_imp.y1_gender_m4];

    
    Sims_Out.P_sim = [Sims_Out.P_sim, P_rep];
     
end
 
 


